Shhh! i want to enjoy the movie!
One of the most annoying things one can encounter is being in the cinema, trying to enjoy a movie and having people around that talk. In this study we will try to link talking to the cinema with personality traits. Specifically, we asked participants to report how often they talk to the cinema in a 7 point likert scale and we administered questionnaires regarding agreeableness and consciousness which are part of the big-five perosnality traits. In addition, we administered questionnaires regarding the dark triad which refer to narcissism - Machiavellianism and psychopathy. We mesured also attention impulsivity which is related to psychopathy.
Our main hypotheses are:
Hypothesis 1: Narcissism will be related to talking to the cinema. Hypothesis 2: Psychopathy will be related to talking to cinema via a attention impulsivity Hypothesis 3: Agreeableness and consciousness will be relate to talking to the cinema.
Load libraries, set paths, load data
Libraries and paths
# Libraries
library(psych) # for descriptive statistics and exploratory factor analysis (also needs *GPAroation* installed, will load it automatically, but not give error if it's not installed; so make sure you install this together with dependencies)
library(GPArotation)
library(dplyr) # for data handling
library(data.table) # for re-arranging outputs of factor analysis to make nice plots (could be done in dplyr now, but not easily when originally coded)
library(ggplot2) # for plotting
library(ggpubr) # for plotting histograms
library(PerformanceAnalytics) # for plotting correlations
library(ggcorrplot) # for plotting correlations
library(outliers) # different test to detect outliers
library(lavaan) # for confirmatory factor analysis
library(ltm) # for cronbach a
library(performance) # for inter-item correlation
library(lavaanPlot) # for plotting lavaan objects
library(semTools) # for measurement invariance
library(brms) # for bayesian regression
# Set paths
plotFolder = file.path(basepath,'/results')
statsFolder = file.path(basepath,'/figures')Read file
#read csv file with raw data
setwd("~/Documents/cinema")
df <- read.csv(file="cinemaDat.csv", header=TRUE, sep=";",dec=",")
#Mahalobi’s distance - detecting outliers First check for multivariate outliers using mahalanobi’s distance.
The Mahalanobis distance is the distance between two points in a multivariate space. It’s often used to find outliers in statistical analyses that involve several variables. To determine if any of the distances are statistically significant, we need to calculate their p-values. The p-value for each distance is calculated as the p-value that corresponds to the Chi-Square statistic of the Mahalanobis distance with k-1 degrees of freedom, where k = number of variables. Typically a p-value that is less than .001 is considered to be an outlier.
mahalanobis distance
my.data<-df %>% dplyr::select(ATTENTION, EXTR, AGR, CON, MAC, NARC, PSYCH, AffEmp, CoEmp)
#check cor with outliers
#create new column in data frame to hold Mahalanobis distances
my.data$mahal <- mahalanobis(my.data, colMeans(my.data), cov(my.data))
#create new column in data frame to hold p-value for each Mahalanobis distance
# df= k-1, k= number of variables that you use for testing
my.data$p <- pchisq(my.data$mahal, df=8, lower.tail=FALSE)
#for 8 df the critical value is 26.13
#drop outliers
df<-slice(df, -c(299, 257, 329,10, 131, 47))
In our sample, six participants can be considered as multivariate outliers according to mahalabis’ distance!
Plots and Correlations
Normal Distribution
Next we will proceed by checking if our data follow the normal distribution. We will do this by checking skeweness and kirtosis and by also checking the qq plots for each variables. For sample sizes greater than 300, depend on the absolute values of skewness and kurtosis without considering z-values. Either an absolute skew value larger than 2 or an absolute kurtosis (proper) larger than 7 may be used as reference values for determining substantial non-normality. For more information see: Statistical notes for clinical researchers: assessing normal distribution (2) using skewness and kurtosis (Kim, 2013)
Correlations
Next we will check the correlations. As our depended variable (talking to cinema) is ordinal we will use sperman sho as method. We will also plot the correlation, scatter plots and histograms of the variables of intererest.
Descriptives statistics, plots and correlations
my.data<-df %>% dplyr::select(food, cell, Talk, ATTENTION, EXTR, AGR, CON, MAC, NARC, PSYCH, AffEmp, CoEmp)
describe(my.data)
#plot qq plots for normality
qq.plots<-list()
qq.plots$ATTENTION<-ggqqplot(df$ATTENTION, ylab = "Attentional Impulsiveness")
qq.plots$EXTR<-ggqqplot(df$EXTR, ylab = "Extraversion")
qq.plots$CON<-ggqqplot(df$CON, ylab = "consientousness")
qq.plots$AGR<-ggqqplot(df$AGR, ylab = "Agreeableness")
qq.plots$MAC<-ggqqplot(df$MAC, ylab = "Maciavelism")
qq.plots$NARC<-ggqqplot(df$NARC, ylab = "Narcissism")
qq.plots$PSYCH<-ggqqplot(df$PSYCH, ylab = "Psychopathy")
qq.plots$AffEmp<-ggqqplot(df$AffEmp, ylab = "Affective Empathy")
qq.plots$CoEmp<-ggqqplot(df$CoEmp, ylab = "Cognitive Empathy")
plot.normality<-ggarrange(qq.plots$ATTENTION,qq.plots$EXTR,qq.plots$CON,qq.plots$AGR,qq.plots$MAC,qq.plots$NARC, qq.plots$PSYCH,qq.plots$AffEmp, qq.plots$CoEmp, nrow = 3,ncol=3, hjust = 0)
plot(plot.normality)
#cordata
cordata<-(my.data)
#plot correlation
chart.Correlation(cordata, method = "spearman", histogram = TRUE, pch=19)
cordata<-(my.data)
Hypotheses Testing
In the next step we will check our hypotheses.
Hypothesis 1.
The first hypothesis states that narcissism will be related to talking in cinema. To test this hypothesis we will use Bayesian regression. We will use talking as depended variable and narcissism as the regressor. Moreover, we will control for covoriates of no interest such as age, gender and education. We will use education as monotonic predictor. Since talking is an ordinal variable we will use family = “cumulative”. We will assess if the one-tailed 95% Bayesian credible intervals of the regression weights exclude zero (significant) or include zero (not significant).
# Settings for BRMS
niter=8000
nchains=4
adaptdeltas=0.9
library(brms)
h1.<-brm(formula=Talk
~1 + NARC+gender+age+mo(educ), data=df,save_all_pars = TRUE, family='cumulative',iter=niter,chains=nchains,cores=nchains,control=list(adapt_delta=adaptdeltas))
## Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
## clang -arch arm64 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG -I"/Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/Rcpp/include/" -I"/Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/RcppEigen/include/" -I"/Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/RcppEigen/include/unsupported" -I"/Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/StanHeaders/include/src/" -I"/Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/StanHeaders/include/" -I"/Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/RcppParallel/include/" -I"/Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/rstan/include" -DEIGEN_NO_DEBUG -DBOOST_DISABLE_ASSERTS -DBOOST_PENDING_INTEGER_LOG2_HPP -DSTAN_THREADS -DBOOST_NO_AUTO_PTR -include '/Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp' -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1 -I/opt/R/arm64/include -fPIC -falign-functions=64 -Wall -g -O2 -c foo.c -o foo.o
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/RcppEigen/include/Eigen/Core:88:
## /Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:1: error: unknown type name 'namespace'
## namespace Eigen {
## ^
## /Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:16: error: expected ';' after top level declarator
## namespace Eigen {
## ^
## ;
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/RcppEigen/include/Eigen/Dense:1:
## /Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/RcppEigen/include/Eigen/Core:96:10: fatal error: 'complex' file not found
## #include <complex>
## ^~~~~~~~~
## 3 errors generated.
## make: *** [foo.o] Error 1
summary(h1.)
## Family: cumulative
## Links: mu = logit; disc = identity
## Formula: Talk ~ 1 + NARC + gender + age + mo(educ)
## Data: df (Number of observations: 352)
## Draws: 4 chains, each with iter = 8000; warmup = 4000; thin = 1;
## total post-warmup draws = 16000
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept[1] -0.35 0.82 -1.90 1.30 1.00 10261 9701
## Intercept[2] 1.18 0.82 -0.38 2.83 1.00 10260 10070
## Intercept[3] 1.95 0.82 0.38 3.61 1.00 10251 9510
## Intercept[4] 2.86 0.83 1.28 4.56 1.00 10138 9788
## Intercept[5] 4.01 0.86 2.38 5.75 1.00 9831 8948
## Intercept[6] 5.49 0.95 3.70 7.46 1.00 9788 9120
## NARC 0.35 0.11 0.13 0.57 1.00 14181 10504
## gender 0.27 0.27 -0.26 0.81 1.00 12636 10560
## age -0.03 0.01 -0.05 -0.01 1.00 12310 11912
## moeduc 0.43 0.44 -0.16 1.41 1.00 2395 5430
##
## Simplex Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## moeduc1[1] 0.13 0.13 0.00 0.47 1.00 4425 3085
## moeduc1[2] 0.07 0.10 0.00 0.39 1.00 3929 8358
## moeduc1[3] 0.07 0.10 0.00 0.37 1.00 4220 7892
## moeduc1[4] 0.08 0.10 0.00 0.37 1.00 5373 8389
## moeduc1[5] 0.37 0.23 0.01 0.80 1.00 3336 3031
## moeduc1[6] 0.27 0.20 0.01 0.72 1.00 7221 8130
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## disc 1.00 0.00 1.00 1.00 NA NA NA
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).Hypothesis 2.
Hypothesis 2 states that psychopathy is related to talking via attentions impulsivity. To test this hypothesis we will use psychopathy as a predictor, talking as independent variable and attention implulsivity as mediator. We will treat the afforementioned variables as latent variables extracted by the observed variables (the items of each questionnaire). Since both the items as well as the dependend variable are ordinal we will use DWLS estimator which is suitable for ordinal data, (see Mîndrilă, 2010; Cheng-Hsien Li, 2016)
To assess good-fitting of the model χ2 and its degrees of freedom (df) were used. For χ2 values associated with p>0,5 were considered good-fitting models, although it has to be mentioned that the p value of this test is sensitive to large sample size. In addition, the root mean square error of approximation (RMSEA) with its 90% confidence intervals (CI), the standardized root mean square residuals (SRMR) and the comparative fit index (CFI) were used. For RMSEA and SRMR values values < 0.08 are acceptable. For CFI values > 0.90 were considered as indicators of good fit (Brown, 2006).
Regarding the mediation, we computed the cross product of the two direct paths coefficient to obtain the indirect path coefficient. Statistical significance level was set at α=0.05.
model2<-'
psych=~Psych1+Psych2_r+Psych3+Psych4+Psych5+Psych6+Psych7_r+Psych8+Psych9
Attentio=~Att1+Att2+Att3_r+Att4+Att5_r+Att6+Att7+Att8
Attentio~a*psych
Talk~b*Attentio
Talk~c*psych
ab :=a*b
total :=c+(a*b)
'
sem.fit2 <- sem(model2, data=df ,estimator = "DWLS")
summary(sem.fit2, fit.measures=TRUE, standardized = TRUE,rsq=TRUE)
## lavaan 0.6-11 ended normally after 68 iterations
##
## Estimator DWLS
## Optimization method NLMINB
## Number of model parameters 38
##
## Number of observations 352
##
## Model Test User Model:
##
## Test statistic 200.146
## Degrees of freedom 133
## P-value (Chi-square) 0.000
##
## Model Test Baseline Model:
##
## Test statistic 1252.229
## Degrees of freedom 153
## P-value 0.000
##
## User Model versus Baseline Model:
##
## Comparative Fit Index (CFI) 0.939
## Tucker-Lewis Index (TLI) 0.930
##
## Root Mean Square Error of Approximation:
##
## RMSEA 0.038
## 90 Percent confidence interval - lower 0.027
## 90 Percent confidence interval - upper 0.048
## P-value RMSEA <= 0.05 0.973
##
## Standardized Root Mean Square Residual:
##
## SRMR 0.062
##
## Parameter Estimates:
##
## Standard errors Standard
## Information Expected
## Information saturated (h1) model Unstructured
##
## Latent Variables:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## psych =~
## Psych1 1.000 0.854 0.499
## Psych2_r -0.456 0.085 -5.347 0.000 -0.390 -0.226
## Psych3 0.994 0.120 8.268 0.000 0.849 0.625
## Psych4 0.856 0.104 8.233 0.000 0.731 0.597
## Psych5 1.490 0.173 8.604 0.000 1.273 0.666
## Psych6 1.059 0.129 8.236 0.000 0.905 0.539
## Psych7_r -0.423 0.086 -4.901 0.000 -0.361 -0.199
## Psych8 0.567 0.097 5.861 0.000 0.484 0.267
## Psych9 1.115 0.134 8.348 0.000 0.953 0.650
## Attentio =~
## Att1 1.000 0.215 0.253
## Att2 1.204 0.263 4.574 0.000 0.259 0.269
## Att3_r -1.691 0.328 -5.157 0.000 -0.364 -0.452
## Att4 2.692 0.499 5.389 0.000 0.579 0.618
## Att5_r -1.100 0.242 -4.544 0.000 -0.237 -0.300
## Att6 1.238 0.273 4.533 0.000 0.266 0.302
## Att7 1.931 0.372 5.191 0.000 0.415 0.444
## Att8 3.082 0.567 5.433 0.000 0.663 0.781
##
## Regressions:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## Attentio ~
## psych (a) 0.080 0.017 4.779 0.000 0.319 0.319
## Talk ~
## Attentio (b) 2.072 0.501 4.135 0.000 0.446 0.299
## psych (c) 0.018 0.078 0.226 0.821 0.015 0.010
##
## Variances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .Psych1 2.203 0.260 8.487 0.000 2.203 0.751
## .Psych2_r 2.814 0.198 14.235 0.000 2.814 0.949
## .Psych3 1.126 0.257 4.379 0.000 1.126 0.609
## .Psych4 0.968 0.188 5.148 0.000 0.968 0.644
## .Psych5 2.034 0.308 6.614 0.000 2.034 0.556
## .Psych6 2.003 0.254 7.879 0.000 2.003 0.710
## .Psych7_r 3.176 0.339 9.376 0.000 3.176 0.961
## .Psych8 3.059 0.266 11.510 0.000 3.059 0.929
## .Psych9 1.240 0.254 4.890 0.000 1.240 0.577
## .Att1 0.675 0.058 11.601 0.000 0.675 0.936
## .Att2 0.859 0.055 15.665 0.000 0.859 0.928
## .Att3_r 0.517 0.050 10.242 0.000 0.517 0.796
## .Att4 0.543 0.079 6.846 0.000 0.543 0.618
## .Att5_r 0.565 0.043 13.078 0.000 0.565 0.910
## .Att6 0.707 0.058 12.132 0.000 0.707 0.909
## .Att7 0.702 0.061 11.440 0.000 0.702 0.803
## .Att8 0.281 0.082 3.430 0.001 0.281 0.390
## .Talk 2.013 0.169 11.918 0.000 2.013 0.908
## psych 0.730 0.127 5.736 0.000 1.000 1.000
## .Attentio 0.042 0.014 3.025 0.002 0.898 0.898
##
## R-Square:
## Estimate
## Psych1 0.249
## Psych2_r 0.051
## Psych3 0.391
## Psych4 0.356
## Psych5 0.444
## Psych6 0.290
## Psych7_r 0.039
## Psych8 0.071
## Psych9 0.423
## Att1 0.064
## Att2 0.072
## Att3_r 0.204
## Att4 0.382
## Att5_r 0.090
## Att6 0.091
## Att7 0.197
## Att8 0.610
## Talk 0.092
## Attentio 0.102
##
## Defined Parameters:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## ab 0.167 0.038 4.418 0.000 0.142 0.096
## total 0.184 0.065 2.815 0.005 0.157 0.106
lavaanPlot(model=sem.fit2, coefs = TRUE, stars = "regress", digits = 2, stand = TRUE)Hypothesis 3
Hypothesis 3 states that Agreeableness and consciousness will be relate to talking to the cinema. To test this hypothesis we will do model comparison. Models give means of putting scientific hypotheses into practice; as a result, model comparison and hypothesis testing are intertwined (Bruno Nicenboim, Daniel Schad, and Shravan Vasishth, 2022).
The Bayes factor reflects a ratio that provides information about how much more likely the observed data are between two compared models. Bayes factor greater than 1 provides evidence in support of one model over the other. We will use Bayes factor to compare the model described in the above-mentioned hypothesis while controlling for covariates of no interesting, against a null model that involves only the intercept without any predictor and only covariates of no interested (age, gender and education treated as monotonic). For more information see vignettes/bayes_factors.Rmd https://easystats.github.io/bayestestR/articles/bayes_factors.html
Also for more information regarding model comparison see:
- Makowski D, Ben-Shachar MS, Lüdecke D. bayestestR: Describing Effects and their Uncertainty, Existence and Significance within the Bayesian Framework. J Open Source Softw. 2019 Aug 13;4(40):1541.
- Vehtari A, Gelman A, Gabry J. Practical Bayesian model evaluation using leave-one-out cross-validation and WAIC. Stat Comput. 2017 Sep;27(5):1413–32.
#fit the null model
null.h3<-brm(formula=Talk
~1 + +gender+age+mo(educ), data=df,save_all_pars = TRUE, family='cumulative',iter=niter,chains=nchains,cores=nchains,control=list(adapt_delta=adaptdeltas))
## Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
## clang -arch arm64 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG -I"/Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/Rcpp/include/" -I"/Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/RcppEigen/include/" -I"/Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/RcppEigen/include/unsupported" -I"/Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/StanHeaders/include/src/" -I"/Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/StanHeaders/include/" -I"/Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/RcppParallel/include/" -I"/Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/rstan/include" -DEIGEN_NO_DEBUG -DBOOST_DISABLE_ASSERTS -DBOOST_PENDING_INTEGER_LOG2_HPP -DSTAN_THREADS -DBOOST_NO_AUTO_PTR -include '/Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp' -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1 -I/opt/R/arm64/include -fPIC -falign-functions=64 -Wall -g -O2 -c foo.c -o foo.o
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/RcppEigen/include/Eigen/Core:88:
## /Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:1: error: unknown type name 'namespace'
## namespace Eigen {
## ^
## /Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:16: error: expected ';' after top level declarator
## namespace Eigen {
## ^
## ;
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/RcppEigen/include/Eigen/Dense:1:
## /Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/RcppEigen/include/Eigen/Core:96:10: fatal error: 'complex' file not found
## #include <complex>
## ^~~~~~~~~
## 3 errors generated.
## make: *** [foo.o] Error 1
summary(null.h3)
## Family: cumulative
## Links: mu = logit; disc = identity
## Formula: Talk ~ 1 + +gender + age + mo(educ)
## Data: df (Number of observations: 352)
## Draws: 4 chains, each with iter = 8000; warmup = 4000; thin = 1;
## total post-warmup draws = 16000
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept[1] -1.57 0.71 -2.90 -0.12 1.00 6682 4240
## Intercept[2] -0.08 0.70 -1.40 1.38 1.00 6566 4421
## Intercept[3] 0.68 0.70 -0.66 2.13 1.00 6657 4574
## Intercept[4] 1.58 0.71 0.23 3.05 1.00 6214 4403
## Intercept[5] 2.72 0.74 1.33 4.24 1.00 6335 4153
## Intercept[6] 4.19 0.85 2.60 5.94 1.00 6469 4653
## gender 0.22 0.27 -0.31 0.75 1.00 14921 11865
## age -0.03 0.01 -0.05 -0.01 1.00 14143 10820
## moeduc 0.38 0.43 -0.16 1.38 1.00 2621 3770
##
## Simplex Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## moeduc1[1] 0.14 0.13 0.00 0.48 1.00 7141 7745
## moeduc1[2] 0.08 0.11 0.00 0.40 1.00 4194 8314
## moeduc1[3] 0.08 0.10 0.00 0.39 1.00 4490 7935
## moeduc1[4] 0.09 0.10 0.00 0.39 1.00 6240 7770
## moeduc1[5] 0.34 0.22 0.01 0.77 1.00 4797 7732
## moeduc1[6] 0.27 0.21 0.01 0.73 1.00 6570 7433
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## disc 1.00 0.00 1.00 1.00 NA NA NA
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
#fit the h3 model
h3.<-brm(formula=Talk
~1 + AGR+CON+gender+age+mo(educ), data=df,save_all_pars = TRUE, family='cumulative',iter=niter,chains=nchains,cores=nchains,control=list(adapt_delta=adaptdeltas))
summary(h3.)
## Family: cumulative
## Links: mu = logit; disc = identity
## Formula: Talk ~ 1 + AGR + CON + gender + age + mo(educ)
## Data: df (Number of observations: 352)
## Draws: 4 chains, each with iter = 8000; warmup = 4000; thin = 1;
## total post-warmup draws = 16000
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept[1] -3.92 1.00 -5.88 -1.93 1.00 9213 5060
## Intercept[2] -2.39 0.99 -4.32 -0.41 1.00 9706 5261
## Intercept[3] -1.61 0.98 -3.53 0.34 1.00 10070 5293
## Intercept[4] -0.69 0.98 -2.59 1.27 1.00 9621 5302
## Intercept[5] 0.46 1.00 -1.45 2.47 1.00 9326 5109
## Intercept[6] 1.94 1.09 -0.13 4.18 1.00 7751 3005
## AGR -0.21 0.12 -0.45 0.03 1.00 15085 10327
## CON -0.30 0.11 -0.51 -0.09 1.00 12678 10741
## gender 0.20 0.27 -0.33 0.71 1.00 14251 10518
## age -0.02 0.01 -0.04 -0.00 1.00 14470 10257
## moeduc 0.42 0.43 -0.10 1.43 1.00 3079 5146
##
## Simplex Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## moeduc1[1] 0.14 0.13 0.00 0.47 1.00 6253 6638
## moeduc1[2] 0.07 0.10 0.00 0.36 1.00 5758 8898
## moeduc1[3] 0.07 0.09 0.00 0.34 1.00 6220 9164
## moeduc1[4] 0.09 0.10 0.00 0.38 1.00 7207 7549
## moeduc1[5] 0.34 0.21 0.02 0.76 1.00 6932 8011
## moeduc1[6] 0.29 0.21 0.01 0.75 1.00 6187 4948
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## disc 1.00 0.00 1.00 1.00 NA NA NA
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
# Compute Bayes factor
bayes_factor(h3., null.h3)
## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
## Iteration: 5
## Iteration: 6
## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
## Iteration: 5
## Iteration: 6
## Iteration: 7
## Estimated Bayes factor in favor of h3. over null.h3: 38.96313